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ABSTRACT 

We study the probability distribution function (PDF) of the mass density in simulations of super- 
sonic turbulence with properties appropriate for molecular clouds. For this study we use Athena, a 
new higher-order Godunov code. We find there are surprisingly similar relationships between the mean 
of the time-averaged PDF and the turbulent Mach number for driven hydrodynamic and strong-field 
MHD turbulence. There is, however, a large scatter about these relations, indicating a high level of 
temporal and spatial variability in the PDF. Thus, the PDF of the mass density is unlikely to be a 
good measure of magnetic field strength. We also find that the PDF of decaying MHD turbulence 
deviates from the mean-Mach relation found in the driven case. This implies that the instantaneous 
Mach number alone is not enough to determine the statistical properties of turbulence that is out of 
equilibrium. The scatter about the mean-Mach relation for driven turbulence, along with the large 
departure of decaying turbulence PDFs from those of driven turbulence, may illuminate one factor 
contributing to the large observed cloud-to-cloud variation in the star formation rate per solar mass. 
Subject headings: ISM: clouds — ISM: magnetic fields — isothermal — simulations — stars: formation 
— turbulence 



1. INTRODUCTION 

The mechanism that determines the star formation 
rate (SFR) within a molecular cloud (MC) is not well 
understood. Observations show that the SFR per solar 
mass, as measured by the ratio of CO to IR luminos- 
ity, varies by as much as 2-3 orders of magnitude from 
region to region; instead Gao & Solomon (2004) found 
that HCN emission, which measures molecular gas at 
much higher density (>3x 10 4 cm~ 3 ), is the better indi- 
cator of star formation. This suggests that the amount 
of mass to reach high density is the key factor in deter- 
mining the SFR per solar mass within a cloud (McKee 
& Ostriker 2007). Observed non-thermal line widths in 
MCs (Falgarone & Philips 1990) indicate that supersonic 
turbulence may be responsible for creating the high den- 
sity contrasts that lead to clump formation. 

Although magnetic fields have been shown not to 
lengthen appreciably the decay timescale of supersonic 
turbulence (Stone, Ostriker, & Gammie 1998, hereafter 
S98; Mac Low 1999), they do create anisotropy in the 
structures within the medium (Vestuto et al. 2003). 
Observations using Zeeman splitting, such as those de- 
scribed by Crutcher (1999), have found magnetic fields 
in some clouds strong enough (e.g. (3 » 0.04) that they 
cannot be neglected. The magnetic field within a MC is 
typically difficult to detect, motivating the determination 
of new diagnostics (e.g. Heyer et al. 2008). 

The probability distribution function (PDF) of the 
logarithm of mass density can be used to quantify the 
amount of material within a turbulent medium that has 
a given density. As the density is likely to have a strong 
impact on star formation, many groups have investigated 
the properties of such PDFs. Padoan et al. (1997, 
hereafter P97), for example, conducted simulations of 
driven hydrodynamic turbulence, determining a relation 



between the mean of the PDF and the Mach number. 
Most of the mass in the simulated clouds was found to be 
in only a small fraction of the volume, with the width of 
the approximately Gaussian PDF increasing with Mach 
number in a predictable way. Ostriker, Stone, & Gammie 
(2001) investigated the effect of magnetic field strength 
on the mean-Mach relation for decaying turbulence, find- 
ing that the fast magnetosonic Mach number can be used 
to predict only a lower limit on the width of the distri- 
bution. 

Recently, we have undertaken a comprehensive study 
of the properties of supersonic turbulence (Lemaster 
& Stone 2008, hereafter Paper I) with Athena, a new 
higher-order Godunov method (Stone et al. 2008). This 
study represents one of the first applications of Godunov 
methods to the study of supersonic MHD turbulence, 
and therefore represents an important test of previous 
results. In this letter, we investigate the effects of mag- 
netic field strength and Mach number on the PDF. We 
survey a much larger range of Mach number and field 
strength, and at a higher numerical resolution, than has 
been used in previous studies. This provides much bet- 
ter statistics and allows us to constrain the form of the 
relations more tightly. We describe the numerical meth- 
ods used to conduct our simulations in explain the 
focus of our analysis in and present the results in fgj 
Finally, we summarize in <J5] 

2. NUMERICAL METHODS 

The simulations we present here were conducted at a 
resolution of 512 3 with Athena (Gardiner & Stone 2005, 
2008; Stone et al. 2008; Stone & Gardiner 2008), a new 
higher-order Godunov code that exactly conserves mass, 
momentum, and magnetic flux. We solve the equations 
of ideal isothermal MHD, 



g + V.(pv) = 0, 



(1) 
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-^+V-( / ovv-BB + P + B 2 /2) = 0, (2) 
at 

and 

— = Vx(vxB), (3) 

where c s — 1 and P = c 2 p are the isothermal sound speed 
and pressure, on a three-dimensional periodic Cartesian 
grid of length L = 1. We use an approximate nonlin- 
ear Riemann solver (HLLD; Miyoshi & Kusano 2005) for 
our MHD runs and an exact nonlinear Riemann solver 
for our hydrodynamic runs. We integrate our simula- 
tions for roughly 4 dynamical times beyond the satura- 
tion time using a directionally-unsplit van Leer scheme 
(Stone & Gardiner 2008). Here we have defined the dy- 
namical time to be t^ yn = L/(2Ai), where M. = <r v /c s 
is the Mach number and a v = {v^p/p) 1 / 2 is the velocity 
dispersion of the gas, calculated using a mass-weighted 
average. For further details on our numerical methods 
as applied to turbulence, see Paper I. 

We drive our turbulence here in a manner very simi- 
lar to S98. We initialize a uniform, stationary ambient 
medium with density p = 1 and magnetic field parallel 
to the x-axis whose amplitude Bq is fixed by the value 
of (3 = 2c 2 p/B' 2 . We apply divergence- free velocity per- 
turbations before every time step, following a Gaussian 
random distribution peaked at k^Ljli: = 2. Before ap- 
plying the perturbations, we shift them such that no net 
momentum will be added to the grid. We also normalize 
them to give the desired energy injection rate, E/pL 2 c^. 
For our decaying runs, we begin with a snapshot of fully- 
developed turbulence from a driven run and allow it to 
evolve without further energy injection. 

We investigate strong-field MHD (/3 = 0.02) as well as 
hydrodynamic (/3 = oo) turbulence. The magnetic fields 
we use in our simulations correspond to physical val- 
ues of B = 2.0/iG/3" 1 / 2 (T/10K) 1 / 2 (n ff2 /10 2 cm" 3 ) 1 / 2 , 
where T is the temperature and nn 2 is the number den- 
sity of molecular hydrogen. Our simulations are scale- 
free, allowing them to be scaled to any set of physical 
parameters using appropriate choices of p, c s , and L. 
Utilizing the same values given in S98, i.e. L = 2pc, 
n H 2 — 10 3 cm -3 , and T = 10 K, yields an energy in- 
jection rate of E = 0.2L Q with magnetic field strength 
B = 44 fiG. 

3. PROBABILITY DISTRIBUTION FUNCTIONS 

Turbulence in molecular clouds causes converging flows 
where the gas can be compressed to very high densities. 
The PDF of the density tells us the fraction of the mass or 
volume within a cloud that obtains a given density. Since 
self-gravitating clumps can form in the high-density re- 
gions, understanding PDFs is critical for gaining insight 
into the stellar IFM and SFR. If compression and rar- 
efaction events in the turbulent gas within a molecular 
cloud are spatially and temporally independent, the PDF 
of density will have a log-normal distribution (Passot & 
Vazquez-Semadeni 1998). The PDF of the logarithm of 
density, then, will have a normal distribution given by 

f(y)dy = -=L= exp dy 7 (4) 

where y = ln(p/ p), p. is the mean of the distribution, and 
a 2 is the dispersion, with \p\ = a 2 /2. 



Our goal in this Letter is to analyze the relationship 
between the mean of the distribution, which represents 
the median density within the cloud, and the Mach num- 
ber. We have investigated this mean-Mach relation over 
a range of Mach numbers 1.2 < M. < 7.0 for both driven 
hydrodynamic and strong-field MHD turbulence. We 
found that a resolution of 512 3 gave smoother-looking 
PDFs that could be fitted more accurately than those 
from 256 3 simulations, justifying the computational ex- 
pense. Higher resolution also allows us to study scatter 
in the PDF in sub-volumes of the domain. We show in 
Paper I that the simulations have converged by this reso- 
lution; thus, increasing it further is unlikely to affect the 
results. To minimize the influence of intermittency, we 
time-average the PDFs obtained from seven snapshots in 
the saturated state spanning almost 3 dynamical times 
before fitting them. 

Since the tails of the PDFs will deviate from normal 
form due to the effects of intermittency, we fit only bins 
with values of at least 10% of the peak value. We per- 
form a Levenberg-Marquardt least-squares fit with uni- 
form weighting. Once we have obtained the mean, p, of 
the best-fit distribution, we plot it against a function of 
Mach number, £(M) = ln[l + <xM 2 ]. With the appropri- 
ate choice of a, we can obtain a linear relation between 
the PDF mean and this function £(A4). We note that 
Kowal et al. (2007) have recently shown that higher or- 
der statistics can also provide insight into the properties 
of supersonic turbulence, however in this Letter we will 
focus only on the density PDF. 

Observations show a wide range of star formation rates 
for different clouds. To infer from our simulations the 
level of cloud-to-cloud variation we would likely observe, 
we also investigate the mean-Mach relation for regions 
of size comparable to the driving scale. To do this, we 
divide our computational domain into eight equal sub- 
domains, each of resolution 256 3 . We compute the PDF 
in each of these sub-boxes (which we will refer to as 
sub-PDFs) individually and plot their means against the 
Mach number within that sub-box. We do not time- 
average our results, yielding 56 mean-Mach pairs for each 
run. Although the snapshots are at intervals of just un- 
der half a dynamical time on the global scale, the inter- 
val between snapshots is closer to a flow crossing time on 
the scale of the sub-boxes, making the snapshots suffi- 
ciently uncorrelated for our analysis. Substantial scatter 
in the values within a run might help explain the ob- 
served cloud-to-cloud variation in the star formation rate 
in molecular clouds. 

4. RESULTS 

The time-averaged PDFs over the full domain are very 
smooth and approximate Gaussians, particularly in the 
hydro case. Although we have not plotted such a PDF in 
this Letter, a similar example can be seen in Figure 3 of 
Kritsuk et al. (2007). Presented in Table Q] are our PDF 
statistics. The first two columns, respectively, specify 
the magnetic f3 and energy injection rate. The final six 
columns are determined from the sub-volumes and de- 
scribe the error bars in the figures to follow. The values 
oy.s and um,s i n the table are the standard deviation of 
the means determined from the PDFs, not the width of 
the PDFs themselves. 
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TABLE 1 
Driven Turbulence at 512 3 



P 


L 1 
Ej 






Hv 


&V,S 


Mm 


°\M,S 


OO 


500 


6.8 


0.5 


-1.1 


0.13 


0.94 


0.09 


oo 


187.5 


5.0 


0.4 


-0.87 


0.12 


0.76 


0.10 


OO 


70 


3.6 


0.3 


-0.65 


0.08 


0.57 


0.07 


oo 


20 


2.4 


0.1 


-0.40 


0.05 


0.35 


0.05 


oo 


1.75 


1.2 


0.05 


-0.10 


0.02 


0.09 


0.01 


0.02 


500 


6.7 


0.5 


-0.89 


0.12 


0.84 


0.11 


0.02 


187.5 


4.9 


0.3 


-0.71 


0.12 


0.66 


0.11 


0.02 


70 


3.6 


0.2 


-0.56 


0.09 


0.50 


0.09 


0.02 


20 


2.5 


0.1 


-0.34 


0.05 


0.30 


0.05 


0.02 


1.75 


1.2 


0.06 


-0.12 


0.03 


0.11 


0.03 



a ay s an d a M S are the standard deviation of /iy and re- 
spectively. 

4.1. Driven Hydrodynamic Turbulence 

Figure [1] shows the mean-Mach relation found from 
time-averaged PDFs over the full domain for driven hy- 
dro and driven MHD turbulence. We find that a value 
of a — 0.5 in the function £(A4) gives the best linear 
relations for the hydro case. For the volume fraction we 
find 

p v = -0.36 ln[l + 0.5A4 2 ] + 0.10, (5) 
while for the mass fraction we find 



Hm = 0.32 ln[l + Q.5M 2 } -0.10. 



(6) 



The rms residuals for these fits are 8.9 x 10 -3 and 
6.5 x 10~ 3 , respectively. Because the density fluctua- 
tions in subsonic turbulence are not produced by shocks, 
we have no reason to expect these relations to approach 
zero with Mach number. The mean-Mach pairs from the 
time-averaged PDFs fall very close to these relations. 
Over the full range of Mach numbers tested, however, 
the time-averaged means are smaller than those found 
by P97 (nv tM = T0.51n[l + 0.257W 2 ]). To determine 
the magnitude of the effect that the driving may have 
had on the relations, we also compare values determined 
from hydrodynamic turbulence with k p kL/2Tr = 4 (not 
shown), finding that these points fall very close to the 
fc p ki/27T = 2 relations as well. 

Also shown in the figure are the instantaneous mean- 
Mach pairs found from each of the eight sub-boxes. The 
scatter in these points, illustrated by la error bars on the 
plot, is quite significant. The ensemble average of the 56 
points for each run falls at a slightly lower Mach num- 
ber than the value found from the time-average over the 
full domain since the Mach number typically decreases 
on smaller scales. These averages differ slightly from 
those computed from time-averaged full-domain PDFs, 
although they still fall close to the relations found. The 
la scatter in the Mach number is 4%-8%, while in the 
sub-PDF means it is 16%-17% for the lowest Mach num- 
ber run and 10%— 14% for the remaining runs. This scat- 
ter puts some of the instantaneous sub-PDF values in 
the vicinity of the P97 relation. 

The scatter in py and Pm can be converted to scatter 
in the median density within the medium. For the ratio 
of the median to mean density within a cloud, p, we 
can define p + and p_ for clouds where p is la above or 
below the mean. The ratios pv,+/pv,- — ex p(2cr\/ i s) and 
Pm.+ / Pm.- = cxp(2(TA/,s) then generally increase with 
Mach number, ranging from 1.0 to 1.3. 
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ln(l + 0.5 Mach 2 ) 



Fig. 1. — PDF mean versus £(A4) for driven hydro (red 
hexagons) and strong-field MHD (blue squares) turbulence. Only 
time-averaged PDFs over the full domain are used to fit the 
mean-Mach relations (hydro: red dashed; MHD: blue dashed). 
Also shown are the instantaneous values taken from each of the 
eight driving-scale sub-boxes (hydro: gold triangles; MHD: green 
crosses) and the lc error bars on those points (hydro: red; MHD: 
blue). For comparison, the relation from P97 is also shown (black 
dotted). 

4.2. Driven MHD Turbulence 

Figure [T] also includes the mean-Mach relation for 
driven strong-field MHD turbulence, where we continue 
to use a — 0.5. The mean-Mach pairs from the time- 
averaged PDFs again fall very close to these relations, 

0.5M 2 



p v = -0.29 ln[l 
for the volume fraction and 



0.06 



(7) 



Pm = 0.28 ln[l + 0.5A4 2 ] + 0.07 (8) 

for the mass fraction, still yielding means smaller than 
those found by P97. The rms residuals for these fits are 
1.1 x 10~ 2 and 1.6 x 10~ 2 , respectively. 

The instantaneous mean-Mach pairs found from sub- 
boxes have more scatter with a strong magnetic field than 
they did in the purely hydrodynamic case. The la scat- 
ter in the Mach number is 5%-8%, while for the sub-PDF 
means it is 22%-24% for the lowest Mach number run 
and 13%-18% for the remaining runs. Converting this 
to scatter in the median density, we find that the ratios 
pv.+/pv,- and pm.+ / 'f>M,- range from 1.1 to 1.3. This 
again puts some of the instantaneous sub-PDF values in 
the vicinity of the P97 mean-Mach relation. 

The time-averaged mean-Mach relations found for hy- 
dro and MHD differ, as one should expect due to differ- 
ences in the shock-jump conditions. However, the sub- 
PDF values overlap substantially, making them difficult 
to distinguish observationally. The time-averaged rela- 
tions fall a bit less than la (computed from the sub-PDF 
values) apart. 

4.3. Decaying MHD Turbulence 

As it seems likely that molecular clouds are transient 
entities, it may be more appropriate to study the PDF 
of decaying turbulence. In Figure we show the evolu- 
tion of the PDF in a decaying strong-field MHD turbu- 
lence run, initialized from a snapshot of fully-developed 
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Fig. 2. — PDF mean versus £(.M) for decaying MHD turbulence 
(orange solid) as it evolves over time (right to left), contrasted 
with the time-averaged driven values (blue squares) and relation 
found in N4.2I (blue dashed). The decaying run does not follow the 
relation found for driven MHD, even considering the large scat- 
ter (error bars and green crosses) found in the driven sub-PDFs. 
For comparison, the hydro relation from P97 is also shown (black 
dotted). 

turbulence from our highest Mach number driven run. 
Although this snapshot has a full-domain PDF mean 
roughly la more extreme than the time-averaged driven 
relations, this should not affect the results. At first the 
small change in mean as the Mach number decreases 
causes a shallower slope than that of the driven relation. 
Once the mean begins to change appreciably, however, 
the slope becomes much steeper, crossing the driven re- 
lation at roughly Ai — 4.5. Although the slope shallows 
as low Mach numbers are reached, it does so only as the 
means become very small. 

The evolution of this decaying run does not parallel or 
asymptotically approach the driven relation at the same 
magnetic (3. Contrary to what was found by Passot & 
Vazquez-Semadeni (1998) for one-dimensional polytropic 
gas, then, we find that the instantaneous Mach number 
is insufficient to describe the properties of the turbu- 
lent medium when the gas is not in a statistically steady 
state. This is consistent with the findings of Ostriker et 
al. (2001) for the PDFs of decaying turbulence. If MCs 
contain decaying turbulence, which seems likely (Heitsch 
et al. 2006), the PDF may be "out of equilibrium" , mak- 
ing relations obtained from steady state (driven) turbu- 
lence inapplicable. 
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5. SUMMARY 

For both supersonic hydrodynamic and strong-field 
MHD turbulence, we have found a one-to-one correspon- 
dence between the mean of the time-averaged PDF and 
the Mach number. The mean-Mach pairs used to fit these 
relations have very small residuals; however, the mean of 
the PDF at any given Mach number, in both the hydro 
and MHD cases, is smaller than was found by P97 for the 
purely hydrodynamic case. Although there is substantial 
scatter of the mean-Mach pairs computed from instanta- 
neous sub- volumes, the ensemble average of these values 
still falls close to the time-averaged global relation. The 
scatter puts a small fraction of the instantaneous values 
in the vicinity of the P97 relation. 

Although the relations found for hydro and MHD dif- 
fer, the scatter about the mean relation of the instan- 
taneous sub-PDF values creates substantial overlap be- 
tween the two. Since the MHD relation gives means that 
are typically smaller than the corresponding hydrody- 
namic values by only la, it will be very difficult to dis- 
tinguish between the two observationally. We have also 
found that PDFs in decaying MHD turbulence differ from 
those of driven MHD turbulence at the same magnetic /3. 
It would seem that the instantaneous Mach number alone 
does not adequately describe the statistical state of the 
turbulent gas when not in equilibrium. Again, however, 
there is substantial overlap between the equilibrium and 
non-equilibrium values, preventing this diagnostic from 
being used to distinguish between the two. 

The scatter we have found about the mean-Mach rela- 
tion may help explain the large variation in the observed 
SFR per solar mass in molecular clouds. Since MCs are 
likely to be transient entities, relations found from driven 
turbulence may not even be applicable to real clouds. If 
large departures from the mean-Mach relation are in fact 
linked to the large variation in the SFR per solar mass, 
this may indicate that turbulence in MCs is indeed de- 
caying rather than forced. 
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